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Numerical simulations of the time-dependent Swift-Hohenberg equation are used to test predic- 
tions of Cross [Phys. Rev. A 25:1065-1076 (1982)] that Rayleigh-Benard convection in the form of 
straight rolls or of an array of dislocations may be observed in an annular domain depending on the 
values of inner radius n, outer radius n, reduced Rayleigh number e, and the initial state. As ri is 
decreased for a fixed r^ and for different choices of e and initial states, we find that there are indeed 
ranges of these parameters for which the predictions of Cross are qualitatively correct. However, 
when the radius difference r2 — r\ becomes larger than a few roll diameters, a new pattern is observed 
consisting of stripe domains separated by radially-oriented grain boundaries. The relative stabilities 
of the various patterns are compared by evaluating their Lypunov functional densities. 
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(N . I. INTRODUCTION 

A central question of noncquilibrium physics is pattern formation, why only certain stationary patterns in a nonequi- 
Qh \ librium, homogeneous, and continuous medium are observed for fixed parameters and boundaries. In this paper, we 
investigate a particular aspect of this question, namely how does pattern formation depend on the shape of the lateral 
boundaries. We do so by numerically integrating a two-dimensional model of a three-dimensional convecting fluid — 
the Swift-Hohenberg equation — in an annular domain as a function of the inner and outer radii of the annulus 
and of the reduced Rayleigh number e (which measures how strongly the fluid is driven out of equilibrium). Theory 
predicts that, as parameters are varied, two kinds of patterns may be observed, a defect-free straight roll texture 
or an array of dislocations The results we present below confirm some of these predictions and also predict a 
\Q \ new state of radially oriented grain boundaries if the radial width becomes sufficiently large. Future Rayleigh-Benard 
■ experiments with a large Prandtl number fluid near onset in an annular domain may be able to confirm these results. 
Our calculations were motivated by experiments |^-^] and simulations |2j which show that a rich variety of stationary 
patterns are observed just above the onset of Rayleigh-Benard convection in the form of stripes (locally periodic 
convection rolls) disordered by topological defects such as dislocations, focus singularities, and grain boundaries. 
(Section VIII of the review article by Cross and Hohenberg Q discusses many experimental and theoretical results 
for Rayleigh-Benard convection.) The occurrence of defects can be qualitatively understood as a consequence of 
two competing effects. First, near onset the local wave numbers of the convection rolls are constrained by the Busse 
stability balloon ^ - 10 1 to lie in a narrow range about the critical wave number q c . Second, in the absence of horizontal 
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thermal gradients, convection rolls are observed to evolve so as to approach the lateral boundaries at an approximately 
normal angle 0. For the most commonly used experimental geometries of cylinders or boxes, these two constraints 
conflict with one another since it is not possible for the rolls to be everywhere normal to the sidewalls and to wiggle 
sinusoidally on a single length scale. This conflict can be partially resolved by the formation of topological defects 
1^3 which allow a local change in wave vector. 

An insightful analysis of this conflict between lateral boundary conditions and local wave number was given by 
Cross [||, who observed that sufficiently close to the onset of convection, the pattern formation was governed by a 
Lyapunov functional C that decreases monotonically along any given orbit and that was bounded from below. Such 
a functional C implies that all initial conditions relax asymptotically to some stationary pattern and also provides a 
way to order different stationary patterns by their relative stability || . In the double limit 

e^0+ and e 1/2 T -> oo, (1) 

where e is proportional to the reduced Rayleigh number (R — R c )/R c and where V is the aspect ratio (ratio of the 
largest lateral dimension to the fluid depth), Cross showed that the Lyapunov functional L could be decomposed into 
a sum of three terms, a surface term associated with the lateral boundaries, a bulk term for the fluid far from the 
boundaries, and a defect term which depends on the particular defects present. By comparing these three terms with 
one another for various geometries in the limit Eq. (nl), Cross was able to explain rather subtle phenomena such as 
why a pattern of straight parallel rolls should be observed sufficiently close to onset in a large cylindrical or annular 
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domain, as opposed to a pattern of concentric rolls with the same axisymmetry as the domain. However, a careful 
comparison of Cross's calculations with experiment (J] or with numerical simulation has not been carried out and it 
remains unclear to what level of accuracy the limits in Eq. (Q) need to be satisfied for Cross's asymptotic analysis to 
hold. Experiments and numerical simulations can also help to determine when higher-order terms in the amplitude 
expansion destroy the property of variational dynamics, leading to a sustained oscillatory or chaotic time dependence. 

Of particular interest to us are Cross's predictions about possible patterns in an annular domain. For this geometry, 
a fundamental frustration can arise if the rolls are normal to both the inner and outer boundaries since it is then 
not possible for the local wave number to remain close to the critical value q c as a function of radius in the limit of 
a large radial difference r 2 — r\. Cross predicted two possibilities that are illustrated in Fig. 7 of his 1982 paper ||. 
One was that sufficiently close to onset, the bulk term in the Lyapunov functional would dominate, which would favor 
straight rolls that extend throughout the annular cell. Wave number frustration and defects arc then avoided since 
the rolls are not normal to the sidewalls. As a second possibility, Cross predicted that the surface term would be 
more important further above onset in which case the rolls would orient normal to the inner and outer boundaries. 
An average wave number close to the critical value q c could then be arranged by introducing an array of dislocations 
(an example is shown below in Fig. [l]). To our knowledge, these predictions have not been tested experimentally or 
numerically. 

In this paper, we test these predictions by numerical simulations of the two-dimensional time-dependent Swift- 
Hohenberg equation jlj in an annular geometry as a simple opportunity to explore how lateral boundaries affect 
pattern formation. The Swift-Hohcnbcrg equation is a widely-studied rotationally invariant model that undergoes 
a supercritical bifurcation from a zero field state to a stripe state, corresponding to the supercritical bifurcation of 
a conducting motionless fluid to finite-amplitude three-dimensional convection rolls. Sufficiently close to the onset 
of convection, the long-wavelength slow-timescale dynamics of the Swift-Hohenberg equation obeys the same ampli- 
tude equation as the one obtained from the three-dimensional Boussinesq equations that describe a convecting fluid 
quantitatively |Q. The linear stability boundaries of the Swift-Hohenberg equation are also similar to that of the 
Boussinesq equations near onset [pj UJl- Numerical integrations of the Swift-Hohenberg equation in periodic, rect- 
angular, and cylindrical geometries (2^ ^],^|] have shown that its solutions are often in semi-quantitative agreement 
with experiments sufficiently close to onset. It is then often useful to explore some questions of pattern formation first 
with Eq. (||) and later make more careful studies by experiment or by more expensive simulations of the Boussinesq 
equations. 

Using a new time-integration code that we have developed to integrate the Swift-Hohenberg equation efficiently 
and accurately in large annular domains |14j , we have calculated transient solutions and their asymptotic stationary 
patterns for different initial states and for different values of the parameters n, r^, £■ As we describe below, our results 
confirm parts of Cross's predictions but we also find that his analysis was incomplete in that a new state occurs when 
the radius difference T2 — r\ becomes sufficiently large. This new state consists of patches of stripes separated by 
radially-oriented grain boundaries. 

Using a modified version of the code, we have also integrated a Generalized Swift-Hohenberg (GSH) model ]H| , [l^] 
to explore whether the wave number frustration caused by the annular geometry may lead to a sustained dynamics. 
The GSH model couples a second field, the vertical vorticity potential £(t,x,y), to Eq. (^) in such a way that the 
overall dynamics is no longer variational and sustained time-dependent solutions are possible |p^-[l9tl . Our admittedly 
brief study suggests that a sustained dynamics exceeding 30 or more horizontal diffusion times (r^) seems to occur 
for r > 15 while for T < 15, long-lived transients were observed that decayed toward a time-independent state. 
However, further studies will be needed to determine more carefully whether any of the states are truly chaotic (what 
might appear to be chaotic states in a spatially extended system might instead be long-lived transients whose average 



decay time grows rapidly with the system size |20 ) and to characterize the dependence of the dynamics on aspect 
ratio and other parameters. 

The rest of this paper is organized as follows. In Section ||, we summarize details of the numerical method used to 
integrate the time-depend ent Swift-Hohenberg equation in an annular geometry and to estimate the average Lyapunov 



functional (£) . In Section III , we discuss the results of our simulations and compare them to the asymptotic analysis 
of Cross and to the Busse linear stability balloon for the Swift-Hohenberg equation jll[] . Finally, in Section IV 
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summarize our main results. Further details are available in the undergraduate Physics thesis of the first author, on 
which this paper is based @. 



II. METHODS 



In this section, we discuss some details of how the Swift-Hohenberg equation was integrated numerically in a large 
annular domain. The two-dimensional time-dependent Swift-Hohenberg equation jl],|^, ^ 



2 



9*^= [e-(V 2 + l) 2 ] V-V' 3 , (2) 

determines the time evolution of a real scalar field ip(t, r, 9) in a two-dimensional spatial domain. Here dt = d/dt 
denotes the partial derivative with respect to time, V 2 = d 2 + <9 2 = <9, 2 + r~ 1 d r + r~ 2 dg is the two-dimensional 
Laplacian operator, and the parameter e is the bifurcation parameter, such that the zero state "0 = becomes linearly 
unstable at a critical value e c (T) > 0. (As shown in Fig. 3 of Ref. 0, for rolls parallel to a wall this critical value for 
the boundary conditions Eq. (^) below goes to zero as (2/r) 2 in the limit r — > oo so e c is tiny for T > 10.) The relation 
of Eq. (Q) to the amplitude equation of the Boussinesq equations shows that, near onset (i.e., in the limit e — > e+), 
the field t/> is proportional to the temperature deviation T(t,x,y, zq) — 7ii noar (zo) of the temperature field T from 
its linear conducting profile Tii near , evaluated on the horizontal midplane z — zq of a three-dimensional convection 
cell. Positive and negative values of tp can therefore be interpreted as warmer (rising) and cooler (descending) regions 
respectively of a convecting fluid. 

We solve Eq. (||) in polar coordinates (r, 9) on the annular domain 

ri < r < r 2 , < 9 < 2vr, (3) 

where the inner radius r% and outer radius r 2 are prescribed parameters. Since the scaling of parameters that lead to 
Eq. (J2J) normalizes the critical wave number to have the value q c — 1, the effective fluid depth for the Swift-Hohenberg 
equation is d = (l/2)(2w/q c ) = n and the aspect ratio T of the annular cell is 

r= r -^I±. ( 4) 

For boundary conditions, we require that the field ip and its normal derivative d r ip vanish at each point of the inner 
boundary r = n and at each point of the outer boundary r = r 2 : 

if)(ri,0) = O, and (<9 r ^)(n, 6) = 0, i = 1, 2. (5) 

Although these conditions are consistent with the vanishing of the amplitude field A — as derived to lowest-order in 
an amplitude expansion in e || , they do not correspond to the correct boundary conditions of the amplitude A further 
from onset. It is then perhaps best to consider the conditions Eq. (^) as convenient phenomenological conditions that 
numerous numerical calculations [l6| have shown to yield reasonable solutions when compared with large-aspect- 
ratio convection experiments near onset. 

The numerical integration of Eq. (||) with boundary conditions Eq. (||) produces a discrete numerical solution ij)^ 
at successive equally spaced times U = iAt that approximates the unknown analytical solution ip{t, r, 9) on mesh 
points (rj,9k) defined by 

rj = ri +jAr, 0<j<N r , (6) 
9 k = kA9, 0<k<N e . (7) 

Here Ar = (r 2 — r{)/N r and A6> = 2ir/(Ng + 1) are the radial and angular resolutions, and the discrete solution 
converges to the unknown analytical solution at the mesh points in the limit of vanishing space-time resolution: 

il>M->1>(ti,r J ,0 k ) as Ai,Ar,A0^O. (8) 

We discretized Eqs. (^) and (^|) on the interior and boundary mesh points by approximating all spatial derivatives 
with second-order-accurate finite differences in the quantity ipijk as described by Bj0rstad in Ref. pjj |. For example, to 
approximate the biharmonic operator V 4 to second-order accuracy at a given lattice point (rj, 9k), a linear combination 
of 13 lattice values involving up to second-nearest neighbors was used. 

Since an explicit time-stepping method for Eq. 
Cmax(Ar 4 , (ArA6*) 4 ) where C is some constant 



would impose a severe stability constraint of the form At < 
22), we chose to integrate Eq. (|j) with a simple-to- implement 
operator-splitting method that integrates the nonlinear term explicitly and then the linear terms in Eq. (^) implicitly. 
The absence of spatial derivatives in the cubic nonlinear term then implies that the largest time step allowed by the 
algorithm is no longer bounded by some power of the spatial resolution, although it remains restricted by the overall 
stability of the time-splitting algorithm and by accuracy. By making runs with different temporal resolutions for a 
fixed fine spatial resolution, we found that a time step of At = 0.5 provided a reasonable balance between accuracy and 
efficiency and we used this value for most of the integrations described in the next section. The code was sufficiently 
fast on a workstation using a 667-MHz Alpha processor that we could integrate tens of horizontal thermal diffusion 
times (many multiples of T 2 ) for our largest T domains in just a few days and so find approximately time-independent 
asymptotic patterns to good accuracy. 
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The first step of our algorithm integrates the nonlinear piece N[i/j] = — -0 3 explicitly, using a second-order-accurate 
Adams-Bashforth method: 

r=i>i + ^(3N{H-Nbl>i-i])- (9) 

This advances the known field field values ipi at time tj (we temporarily suppress the spatial indices ij to simplify the 
notation) to intermediate field values that we denote by 0*. The field values at time U + \ are then obtained by solving 
the linear pieces in Eq. (j^) implicitly using a backward-Euler method with initial data ip*, leading to the equation 

V 4 + 2V 2 + l + ^-e) Vi+i = (10) 

with the boundary conditions Eq. (|^) imposed on ipi+i- This constant-coefficient generalized linear biharmonic 
problem can be solved efficiently with a second-order-accurate fast-direct method invented by Bj0rstad [|l],^,|2j] 
which has a nearly optimal computational complexity of 0(N r Ng \og(Ng)). We used the public Fortran77 version of 
the annular Bj0rstad solver available through netlib pq| . 

The input data for the computer code are then the parameters e, n, r"2, N r , Ng, At, the total integration time T, 
and some initial discrete field ipojk on the mesh Eq. (p|) that satisfies the discrete form of the boundary conditions 
Eq. (||). Previous studies Q have shown that a reasonable spatial resolution requires at least six mesh points per half 
roll (distance of n) or Ar < tt/6 w 0.52. For fixed r\ and r2, we therefore chose the spatial resolutions N r and Ng 
to provide or to exceed this resolution in the radial and azimuthal directions. (Although the resolution is uniform in 
the r and 9 coordinates, in an annular cell geometry the field ip is approximated more accurately near the inner radius 
since the density of angular mesh points is higher there.) Many of our runs were tested at higher spatial and temporal 
resolutions (usually by a factor of 2 or 4) and our results were verified not to change with the higher resolution. 

We used four different kinds of initial conditions ipo( r ,6) — ip(t — 0,r,9) to explore the dependence of the final 
pattern on the initial state, 

■00,1 = acos(k e 9), ^0,2 = asm(k r r + 0), , . 

"00,3 — a sin (k x T cos 9 + <f\ , ipo,4 = a7]. 

These describe respectively radially oriented rolls (ipo,i), azimuthally oriented (concentric) rolls (^0,2), straight rolls 
(■0o,3), and random noise (^0,4 )■ (We also occasionally added small-amplitude noise to the roll initial conditions to 
break possible symmetries.) The parameter a denotes the initial amplitude (often chosen to be O.l-^/e), the parameters 
kg, k r , and k x are specified wavenumbers, is some specified phase, and the 77 are uniformly distributed random 
numbers in the interval [—a, a]. Although the second step Eq. ( 10|) of the operating-splitting method automatically 
makes the solution ip(t + At) satisfy the discrete boundary conditions Eq. (g), we found that convergence to a smooth 
solution was enhanced if we forced the initial states Eq. to satisfy the boundary conditions by setting the field tpo 
to zero within a distance Ar = n of the inner and outer radii. 

The code then generated spatiotemporal fields ipijk on the polar jfc-mesh at successive times U = iAt. We visualized 
the fields with contour plots at positive and negative contours of magnitude ±(1/2) max^fc IVfy'fcl- We also plotted 
three scalar time series that proved useful when determining whether some state had become stationary. One time 
series was the local field value at the midpoint of the annulus at angle 9 = 

Sl ( t )=V,(t,^±^,o). (12) 

The second series was a global quantity, the spatially averaged mean-square field (tp 2 ) (t) 

s 2 (t) = j £ dr J 2 Jd9 r m, r, 9)f (13) 

where A = n(r\ — r\) is the area of the annular domain. Near onset, one can show that Eq. ( |T3] ) is linearly related 
to the global vertical heat transport (Nusselt number) across the fluid layer and so a plot of S2 (t) indicates how heat 
transport depends on the spatiotemporal dynamics. As a third time series, we recorded the instantaneous value of 
the spatially averaged Lyapunov functional (£[0]) 

s 8 (*) = »]>(*) =^], (14) 
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where [||: 



Cty] = f "dr f ddr \-eip 2 + ^ + {(V 2 + 1)V>} 2 ] • (15) 
Jr x Jo l ^ J 

Because the mesh is uniform with respect to the r and 9 variables, the integrals in Eqs. ( |l3| ) and ( |l4| ) could be approx- 
imated with a simple two-dimensional version of the rectangle rule for approximating one-dimensional integrals p6| . 
The integrands were first approximated at the mesh points Eq. (^) using second-order accurate finite-difference ap- 
proximations for the spatial derivatives, and then the mesh values were accumulated over all the mesh points. For 
many initial conditions, we verified that the time series S3(t) in Eq. (|lj) decreased monotonically with time provided 
that the space-time resolution was sufficiently fine, as should be the case for the analytical functional associated with 
Eq. (Q) and with boundary conditions Eq. (^|). 

After observing many patterns and their time series, we chose empirically to define a pattern to be stationary if 
each of the two global time series Eq. (|l3|) and Eq. (|lj) had converged to three significant digits. In fact, many of the 
patterns continued to evolve extremely slowly even after many horizontal diffusion times (since the phase of the stripe 
amplitude is diffusing in a large domain), and we saw changes in the fourth significant digit of the quantities (tp ) 
and (£). However, when we integrated by a factor of ten longer for certain representative patterns (nearly 100 T 2 in 
some cases), we found no significant changes in the patterns except for tiny translations or rotations of the overall 
roll structure. Although future calculations of stationary states may benefit by using more sophisticated algorithm 
such as a Newton method [^6| applied to the time-independent form of Eq. (0) or by using an annular version of a 
recently developed fully-implicit method flj|| , the present results should provide a good starting point for comparing 
theory with experiment. 

III. RESULTS AND DISCUSSION 
A. Overview 

In this section we summarize results of numerical integrations in annular domains of the Swift-Hohcnberg equation 
Eq. (^|) for different choices of the bifurcation parameter e, inner and outer radii r\ and r 2 , and initial conditions. 
Our results confirm some predictions of Cross Q but also show reveal new patterns that arise when the aspect ratio 
Eq. (|]) becomes sufficiently large. We also mention a few results for the Generalized Swift-Hohenberg equation which 
permits time-dependent nontransient states. 

The space of parameters e, n, r 2 , and initial data is too large to explore systematically so we let the predictions of 
Cross, as well as practical computational constraints guide us in our choice of parameters. Our interest was first to 
explore the region defined by the limits Eq. (^) for which the amplitude equation should give an accurate description 
of the dynamics. For the annular domain Eq. (^), Eq. (|l|) plus the assumption of Cross that the annular aspect ratio 
is not too large imply the following constraints 

rjv^»l, (r 2 - nK/i > 1, (16) 

ri 

and we would like to satisfy these inequalities as best as possible. A second goal was to test the specific prediction of 
Cros s B that a texture involving roll dislocations becomes more stable (has a lower Lyapunov functional value C in 
Eq. (|l4j)) when the following inequality is first satisfied: 

eV4 > r2-n 
n 

For fixed radii T\ and r 2 , Eq. (117) can be tested by making runs for e above and below the ratio [(r 2 — fi)/fi] 4 for 
the four initial conditions Eq. (Ill ) , followed by comparing the average value of the Lyapunov functional (C) for the 
final stationary states. 

Another constraint on the choice of parameters e, r\, and r 2 for an annular domain comes from the Busse stability 
balloon for the Swift-Hohenberg equation (see Fig. 31 in Ref.Q). For locally parallel rolls that are normal to the inner 
and outer walls of an annulus, one might expect a pattern change when the local wavenumber q(r) as a function of 
radius crosses a linear stability boundary. For example, if the local wave number q near a fixed outer radius r 2 is as 
small as possible consistent with the zig-zag stable region q > 1, then the local wave number q(r\) enters the Eckhaus 



unstable region q > 1 + y/e/12 when 
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n < — r A=- (18) 
' - 1 + ^712 

For e = 0.1 and fixed T2 = 80 (two parameters used in the runs below), the local wavenumber q{r\) becomes Eckhaus 
unstable when r\ < 73.3. The corresponding aspect ratio T w 2.1 is so small that the critical value e c for the instability 
of the zero- field state in a finite system exceeds e = 0.1, i.e., as soon as a finite amplitude radial stripe pattern can 
form, it is already Eckhaus unstable. A larger outer radius r-i would therefore be needed to test the implications of 
Eq.©. 

Given the above constraints and our computational resources, we chose to fix the outer radius at the value T2 = 80 
(about 25 half-rolls) and study patterns with a decreasing sequence of inner radii, from n = 60 to r\ = 10 in steps 
of 10. This corresponds to increasing the aspect ratio T from 6.4 to 22.3 in steps of Ar = 3.2. Eq. ( |l6| ) then 
implies that we can make comparisons with Cross's prediction Eq. ( \n\ ) for T < 10. For the different choices of rj 
for fixed r%, we studied just two values of the bifurcation parameter, e = 0.1 and e = 0.5. Earlier calculations of 
the Swift-Hohenberg equation in rectangular domains S show that barriers to dislocation motion appear for e larger 
than about 0.5 so these two values probe the "near onset" and "not so close to onset" regimes. 



B. Stationary Patterns 

We begin our discussion of results with Fig. [j], which shows two asymptotic stationary patterns near onset (e = 0.1) 
obtained by integrating the Swift-Hohenberg equation for more than 15 horizontal diffusion times (units of t/j = T 2 ) in 
a small aspect ratio T ~ 6.4 annular domain of outer radius = 80. The initial conditions for panels (a) and (b) were 
respectively radially oriented rolls and straight rolls (V'0,1 and "00,3 i n Eq- (|Tl|)). (Random initial conditions lead to the 
same kind of state so a pattern with dislocations seems to be the most accessible basin of attraction.) The asymptotic 
patterns in both cases consist of approximately straight radially oriented rolls disrupted by four to six dislocations. 
The states in (a) and (b) have an identical average Lyapunov functional (£) = 3.40 x 10 -3 to three significant digits 
and so are extremely close in "energy" . Initial conditions consisting of concentric rolls of high symmetry and a little 
noise evolve into nonlinear stationary concentric rolls which therefore constitute another basin of attraction, although 
a small one. The average value of the Lyapunov function (£) = 0.0819 is higher than that for any of the patterns 
with dislocations so the concentric rolls are not preferred. 

The value e = 0.1 in Fig. [j] substantially exceeds the value \{v2 — ri)/ri] 4 » 0.012 predicted by the criterion Eq. ( |l7| ) 
with r2 = 80, r\ = 60 and so Fig. [I] is consistent with Cross's prediction that states with dislocations are preferred in 
that they have a lower average Lyapunov functional. In fact, because the threshold e ~ 0.012 defined by Eq. (|l7]) is 
smaller than the critical value e c for the onset of convection in such a small aspect ratio cell, we could not find a pattern 
of straight rolls for these same parameters to compare Lyapunov functional values. For the larger value of e = 0.5, 
which also satisfies Eq. (|l7]), asymptotic patterns similar to Fig. |l| are observed for initial conditions consisting of 
radially oriented rolls or of random values. These patterns further from onset have an average Lyapunov functional 
and average square value of (£) = —0.0875 and (ip 2 ) = 2.79 x 10~ 5 . But for an initial condition consisting of straight 
rolls at the critical wave number q = 1, Fig. |2| shows that a dislocation- free pattern of approximately straight rolls 
can occur, although the stripes are only approximately straight where they are nearly parallel to the inner and outer 
boundaries. This pattern has average values of (£) = 0.0819 and (■0 2 ) = 2.74 x 10~ 4 and so is not preferred, again in 
accordance with the criterion Eq. ([17|). 

We now increase the aspect ratio from T = 6.4 to 12.8 by decreasing the inner radius from r% = 60 to 40 for 
fixed T2 = 80. Cross's criterion Eq. ( |l7| ) for a texture with dislocations to be preferred is now no longer expected to 
hold since the annulus is many rolls wide and more complex textures can be expected. Fig. || shows two examples 
of the approximately stationary patterns observed after long times near onset (e = 0.1), starting from two different 
initial conditions, radially oriented rolls in (a) and straight rolls in (b). The patterns have some similarities to Fig. |l| 
in that the rolls are largely radially oriented but they now bend at a substantial angle as they traverse the outer to 
inner boundaries and a novel defect structure is observed in the form of four or five radially-oriented grain boundaries. 
These and other runs suggest that there are many different stationary states possible for the specified parameter 
values, all differing slightly in the position and number of the grain boundaries and possibly with the presence of a 
few other defects. For example, in addition to four grain boundaries Fig. ^(b) shows two dislocations near the top 
and bottom of the inner boundary. According to the values of £, Fig. ||(a) is slightly preferred over (b) but the lack 
of symmetry suggests that neither is the global minimum for Eq. (|l5|). 

For this same T — 12.8 cell, we carried out several runs further from onset for the value e = 0.5. Straight roll 
initial states lead to a straight-roll stationary state (not shown) similar to Fig. || with average values of (£) = —0.176 
and (ip 2 ) = 3.26 x 10~ 5 . Radially-oriented and random initial conditions lead to patterns similar to Fig. [l] but with 
many more defects (with values (£) = —0.176 and (ip 2 ) = 3.19 x 10~ 5 ). Although the straight-roll state has a 
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lower average Lyapunov value, it is not accessible from typical initial conditions because of the barriers to dislocation 
motion. 

Fig. [| shows two stationary patterns near onset (e = 0.1) in a still larger T = 19.2 cell, while Fig. |5| demonstrates the 
stationarity of such patterns by showing how the time series for the average Lyapunov functional (£) and mean square 
field (ip 2 ) have converged to constant values. A new feature in Fig. ^(a) is the appearance of a focus singularity for the 
first time (near the upper right side of the annulus). Foci are a common feature in large cylindrical and rectangular 
patterns near onset but evidently don't occur in annular domains until the aspect ratio is sufficiently large. The 
stationary state Fig. ||(b) arises from straight-roll initial conditions. Again there is a new feature with this larger 
aspect ratio in the form of grain boundaries next to the top and bottom outer boundary. Such grain boundaries are 
a common feature in large rectangular cells near onset (for example, see Fig. 9 of Ref. since rolls parallel to a 
sidewall are known to be unstable near onset to transverse rolls. The defected and straight roll patterns in Figs. ||(a) 
and (b) have identical average Lyapunov functionals to three significant digits and so are equally stable. However, 
the straight roll state empirically has a tiny basin of attraction and so is not accessible for most initial states. 

We conclude this subsection with Fig. which shows two representative patterns in the T = 19.2 cell further 
from onset for e = 0.5. Fig. ^|(a) grew out of straight roll initial conditions and is a nearly perfect array of stripes; 
the instability near sidewalls towards transverse rolls does not occur at this value of e. Fig. 0(b) shows a state that 
grew out of small amplitude random initial conditions. The highly disordered lamellar- like pattern is similar to that 
observed for the Swift-Hohenberg equation in a large rectangular cell for e = 0.5 and is a consequence of larger 
barriers to the motion of dislocations. Such frozen disordered states have not been observed experimentally and are 
likely an artifact of the Swift-Hohcnberg equation, which is most physically relevant to Rayleigh-Benard convection 
in the limit e — > and large Prandtl number. 



C. Dynamical States of the Generalized Swift-Hohenberg Equation 

Convection experiments close to onset of room-temperature pressurized gases (Prandtl number 1) have revealed 
many different sustained time-dependent states of which spiral defect chaos (SDC) Jl^.^7] [h] has been perhaps the 
most intriguing and thoroughly studied. Several experimental features of SDC remain poorly understood, e.g., why 
the state occurs in the first place, why its statistical properties are insensitive to the aspect ratio of a cell (above some 
minimum value of T), and why, to the contrary, the critical value egDC above which SDC is first observed (for a fixed 
experimental cell) is sensitive to the size and geometry of the convection cell. 

To determine whether SDC may be modified by an annular geometry and whether the wave number frustation of 
an annular cell may itself be a source of sustained dynamics, we have carried out a few exploratory calculations with 
a Generalized Swift-Hohenberg model |Ll| , |l5| for which prior numerical calculations in large square periodic domains 
have produced spiral-defect-chaos like states |l^,|2^,^2|. We used a model of the form 



d t ip = [e - (V 2 + l) 2 ]^ - ^ 3 - g(V • V)?/>, (19) 

(-77V 2 + c 2 ) V 2 C = z ■ [V (V 2 </>) x Vip] , (20) 

U = Vx(Cz), (21) 

with boundary conditions Eq. (|5|) for the field ip and boundary conditions 

CM = C(ra) = 0, d r (\ ri = d r C\ r2 = 0, (22) 



for the vorticity potential field (,{t,r,Q). Because Eq. (|2C]) is linear in £ and involves again a constant-coefficient 
generalized biharmonic operator, it was straightforward to modify the numerical algorithm of Section |0] to integrate 
Eqs. @- @ in time. 

For all of our runs, we used the parameters of Ref. |^9| 

e = 0.5, 7] = 1, g = 50, c 2 = 2, (23) 

for which long-lived spiral-defect-chaos- like states were observed with periodic boundary conditions. We then fixed 
the outer radius to be ri = 80 and made several runs, each starting with random initial conditions, for values of 
the inner radii from r\ — 60 to r\ — 10 decreasing in steps of —10. (This corresponds to aspect ratios between 6.4 
and 22.3 in steps of Ar « —3.2.) For T < 13, we found that all states eventually became stationary (using the time 
series of (ip 2 ) as the criterion) with the state T = 12.7 taking the long time of 43 77, to become time independent. 
For r > 13, the states were still dynamic up to our longest integration times of 50 although there was evidence 
both visually and from time series that the states were coarsening and slowly evolving to a state other than SDC. 



7 



A representative example of the spatiotemporal dynamics of the GSH model Eqs. (|19|)- ( |2l| ) in a large aspect ratio 
r = 22.3 annular cell is shown in Fig. fj]. The contour plots show the spatial structure of the field -0 at times t = 4, 12, 
18, and 24 Th- At earlier times (panel a), there are a few spirals but these eventually annihilate leaving a pattern that 
is more slowly evolving, and with nearly half the pattern consisting of foci. Two associated time series are shown in 
Fig. ^. They suggest that the dynamics is more complex for about 10 horizontal diffusion times, then becomes slower 
and simpler, unlike the experimentally observed SDC state. The approximately periodic behavior of the midpoint 
field value in Fig. |8fo for t > f 5 th is likely caused by a slow translation of stripes past the observation point and does 
not correspond to a true time-periodic behavior. These time series suggest that the pattern may eventually become 
time independent. 

Because the approximations relating the GSH model to the Boussinesq equations are less well justified than those 
that relate the Swift-Hohncberg model near onset, especially concerning sustained chaotic states, it would be inter- 
esting to continue these calculations not with further integrations of the GSH model but with integrations of the 
three-dimensional Boussinesq equations in an annular domain, and to test the various calculations with experiments 
in large-aspect-ratio annular domains. 

IV. CONCLUSIONS 

The previous sections may be summarized as follows. We have developed and applied a new computer code to 
study pattern formation formation near onset of the Swift-Hohenberg model in annular domains of varying aspect 
ratios. As pointed out by Cross ||, an annular geometry is interesting because there is an inherent conflict near onset 
between the tendency of rolls to be normal to the sidewalls and the tendency for the rolls to be straight in the bulk. 
We have investigated Cross's predictions for annular cells of varying aspect ratio 6.4 < V < 22.3, near onset with 
e = O.f , and further from onset with e = 0.5. 

For an annular domain of rather small aspect ratio (r < 5), the predictions of Cross hold and stationary states with 
dislocations are preferred (have a lower Lyapunov functional density) compared to straight rolls. As the aspect ratio T 
becomes larger, new patterns are observed that are characterized by radially oriented grain boundaries. In nearly 
all cases, for a given geometry and reduced Rayleigh number e, many different stationary states are possible, usually 
differing only slightly in the value of the Lyapunov density or mean square field (Nusselt number). The pattern that 
produces the global minimum of C for given parameters values e, r%, and r 2 is not known but we conjecture to be 
a highly symmetric arrangement of defects. Straight-roll patterns are observed only when starting from straight-roll 
initial conditions. Such patterns can have a lower Lyapunov density than the patterns with radially oriented rolls but 
evidently have such a small basin of attraction that they are not observed starting from most initial conditions. 

We also explored the possibility of sustained time-dependent states using a Generalized Swift-Hohenberg model, for 
parameters such that a long-lived spiral defect state is observed in a large aspect ratio periodic square |^] . However our 
admittedly incomplete study revealed only long-lived complex spatiotemporal transients. Further analysis, preferably 
by experiment or by integrations of the Boussinesq equations, will be needed to determine the effect of an annular 
domain on time-dependent states. 

We hope that the above results will stimulate experiments to test the specific calculations reported here, especially in 
large Prandtl number fluids near onset for which the Swift-Hohenberg equation should provide a reasonable description. 
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FIG. 1. Stationary numerical solutions of the Swift-Hohenberg equation Eq. ^ with boundary conditions Eq. (JE]) in an 
annular domain of outer radius V2 = 80, inner radius r± — 60, aspect ratio F ~ 6.4, and bifurcation parameter e = 0.1. Part (a): 
A pattern at time t = 17 th (units of horizontal diffusion times th) that grew out of the initial state i^o = 0.3 cos(100#). Part (b): 
a pattern at time 19 Th that grew out of a straight-roll initial state ipo = 0.3 cos(r cos(#)). For both patterns, {£.) [tp] = 3.40 x 10~ 3 
and {ip 2 ) = 4.43 x 10" 6 . The spatiotemporal resolutions were N r = 66, N g = 1024, and At = 0.5. 

FIG. 2. Stationary pattern at time t = 25 tt, starting from straight roll initial conditions tpo — 0.3 cos(r cos(#)) for the same 
physical and numerical parameters as in Fig. hi but with a larger value e = 0.5. Now (£} = 0.0819 and (4> 2 ) = 2.74 x 10~ 5 . 
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FIG. 3. Two stationary patterns in a larger annular domain for parameters e = 0.1, ri = 40, = 80, T « 12.7, At = 0.5, 
N r = 130 and Ne = 1024. (a): An approximately stationary pattern obtained at time t = 19 th starting from azimuthally 
oriented rolls, with values {£} = —4.50 x 10 3 and {ip 2 } = 3.32 x 10~ 6 . (b): An approximately stationary pattern obtained at 
time t = 22 starting from straight roll initial conditions, with (£) — —4.45 x 10 -3 and {ip 2 ) = 3.38 x 10 -6 . Both patterns 
show new defects consisting of radially oriented grain boundaries. 

FIG. 4. Two stationary patterns in an annular domain for parameters e = 0.1, ri = 20, T2 = 80, V ~ 19.2, At — 0.5, 
N r = 258 and N$ = 512. (a): An approximately stationary pattern obtained at time t = 22 starting from radially-oriented 
rolls, (b): An approximately stationary pattern obtained at time t = 22 tk starting from straight roll initial conditions. These 
and similar patterns obtained starting from random initial conditions all have identical average Lyapunov functionals and 
square values to three significant digits, with respectively (£) = —4.86 x 10 -3 and {tp 2 ) = 2.88 x 10 -6 . 

FIG. 5. Time series of the spatially averaged Lypuanov functional (£) (t) (Eq. (II 41)) and mean-square value (ip 2 )(t) (Eq. @) 
starting from low-amplitude random initial conditions for the parameters of FigT~EI. The states are accurately stationary for 
t > 19 th, although it is rather interesting that the two time series do not become stationary at the same time. 

FIG. 6. Two stationary patterns in an annular domain for parameters e = 0.5, r\ — 20, T2 = 80, V ~ 19.2, At — 0.5, 
N r ~ 258 and Ne — 512. (a): An approximately stationary pattern obtained at time t = 2tj, starting from straight rolls, 
with values £ = —0.110 and {ip 2 } = 1-59 x 10 -5 . An essentially perfect straight roll state is attained, (b): An approximately 
stationary pattern obtained at time t = 8 Th starting from small-amplitude random initial conditions. Although there are many 
defects, the rolls are still oriented normal to the boundaries. 

FIG. 7. Contour plots of the field ip(t, r, 9) at times t = 6, 12, 18, and 24 th from an integration of the General- 
ized-Swift-Hohenberg equations Eqs. (|l9|)- ( pi] ) with boundary conditions Eqs. (||) and ( |22[ ) in an annular domain with = 80, 
r\ = 10, r » 22.3 and starting from small-amplitude random initial conditions. The parameters used were e = 0.5, r\ = 1, 
c 2 = 2, g = 50, N r = 514, Ne = 1024, and At — 0.125. The state is time dependent out to 27 th, the longest time observed. 

FIG. 8. Time series of (a), the mean square field (t/i 2 ) and of (b), the midpoint value Eq. ( p^ ) for the time dependent 
solution of the GSH equations shown in Fig. 0. The dynamics is becoming slower and simpler over time. 
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Figure 7. 
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